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ABSTRACT 



Magnetic flux is continuously being removed and replenished on the solar surface. To understand the removal process 
we carried out 3D radiative MHD simulations of the evolution of patches of photospheric magnetic field with equal 
amounts of positive and negative flux. We find that the flux is removed at a rate corresponding to an effective turbulent 
diffusivity, ?7eff, of 100-340 km 2 s _1 , depending on the boundary conditions. For average unsigned flux densities above 
about 70 Gauss, the percentage of surface magnetic energy coming from different field strengths is almost invariant. The 
overall process is then one where magnetic elements are advected by the horizontal granular motions and occasionally 
come into contact with opposite-polarity elements. These reconnect above the photosphere on a comparatively short 
time scale after which the U loops produced rapidly escape through the upper surface while the downward retraction 
of inverse-?/ loops is significantly slower, because of the higher inertia and lower plasma beta in the deeper layers. 
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continuously removed and replenished. The small scale 
fields involved are probably a mixture of flux which 
is left over from the decay of active regions and flux 
generate d in the near-phqtosph e ric layers by local-dynamo 
action (Hagena ar et al.l 120031: IVogler fc Schiisslerl 120071 : 
iPietarila Graham et al.l I2QQ9L 120101 ). This field has been 
described as 'salt and pepper', reflecting small spatial 
scales on which the polarities are mixed. Reconnection 
events take place when the surface magnetic fields of 
opposite polarities, which are continually advected by 
the granular motions, come into close proximity. If we 
consider only a layer near the surface, then a field line 
in the layer can be classified as a U loop if both ends 
of the field line pass through the upper boundary of the 
layer, and as a inverse-// loop if both ends pass through 
the lower boundary. Reconnection of field lines initially 
passing from the bottom of the layer to the top with field 
lines going from the top to the bottom produces pairs of 
U and inverse- U loops. Magnetic tension then causes the 
loops to retract: inverse-// loops downwards through the 
photosphere and U loops upwards. Observations indicate 
that the reconnection occurs above the optical surface, 
so that the removal of flux from the surface invol ves the 
subsequent retraction of U loops ([Kubo et al.ll2010l ). 

In this paper we study the evolution of a small-scale 
mixed polarity field using numerical simulations. We con- 
centrate not so much on individual reconnection events but 



on the statistical properties of the flux removal such as 
how fast does the removal occur, what are the processes 
which determine the rate, what are the dynamics of U and 
inverse-// loops. The observational sign ature of individu al 
reconnection events has been studied by iDanilovid ([.2009) . 

This paper is organized as follows. In Section 2 we 
briefly sketch the numerical code we have used and the 
simulation setup. Section 3 presents the simulation results 
and their dependence on the assumed initial and boundary 
conditions. Section 4 discusses the results, and Section 5 
gives our conclusions. 



2. Simulations 

Our simulations consider the near-surface layers of the Sun 
extending from about 800 km below the optical surface 
to 600 km above. The physics is described by the conti- 
nuity and momentum equations, the induction equation, 
an energy equation which includes radiative transfer, and 
an equation of state which allows for the effe cts of partial 
i oniza tion. The equations are fully set out in I Vogler et al.l 
([2005f h where also the details of the MURaM code we have 
used for this study are described. 

The code has been used to investigate numer- 
ous aspects of solar-surface magnetic activ ity in- 
cluding quiet-S un ma gne tism ([K eller e t al.l 
Kho menko et al.l 120051 : | Vogler fc Schiis slerl 
Pietarila Gr aham et all I2009L 120101: IDanilovic et~aD 



2004 
2007 



2010c. b) pore s (ICameron et al.l I2007D, emerg ing flux 



(ICheung et all 120081 : lYelles Chaouche et all 120091: 



Cheung et al. 2010) , and sunspots ( Schiissler fc Vog ler 



l200a iRempel et al.ll2009bl laf). 

For this study we first consider a reference case with a 
box size of 6 Mmx6 Mm in the x and y (horizontal) direc- 
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tions and 1.4 Mm in the z (vertical) direction with a grid 
resolution of 20.8 km in each of the horizontal directions 
and 14 km in the vertical. The bottom bou ndary condi- 
tion a llows for in- and outflows as described in lVogler et al.l 
(2005). The magnetic field is vertical at the top boundary, 
where the vertical component of the velocity vanishes and a 
free-slip condition (vertical derivatives vanish) is condition 
is imposed on the horizontal velocity components. The box 
is periodic in both horizontal directions. The initial condi- 
tion for our reference case is a height-independent checker- 
board (2 x 2) vertical magnetic field, i.e. 



B z 



-200G 
-200G 
+200G 
1+200G 



for < x < 3 Mm and < y < 3 Mm 
for 3 < x < 6 Mm and 3 < y < 6 Mm 
for 3 < x < 6 Mm < y < 3 Mm 
for < x < 3 Mm 3 < y < 6 Mm. 



B x = Bn, 



o, 



(1) 



superposed on a well developed non-magnetic convection 
simulation. This initial condition is not particularly realistic 
in that we do not expect purely vertical fields to occur over 
a region this large in the quiet Sun. The initial condition is 
instead designed to investigate the manner in which flux in 
the quiet Sun is dispersed, reconnects and is removed from 
the surface. 

We also consider several variations of this setup. We 
study the effect of the initial condition by starting from a 
magnetic field configuration consisting of two stripes (2x1) 
of opposite polarity: 



B z 



J-200G for < x < 3 Mm 
\+200G for 3 < x < 6 Mm. 



(2) 



We also investigate the sensitivity to the upper bound- 
ary condition by comparing the reference case to ones where 
the magnetic field was matched to a potential field at the 
upper boundary, to cases where the top boundary is raised 
by additional 280 km above the optical surface, and to cases 
where the top boundary is open to flows. These closed and 
open upper boundary conditions either completely or par- 
tially reflect outgoing waves. The influence of the reflected 
waves, and of waves otherwise excited in the chromosphere, 
on the retraction or escape of the loops is not studied be- 
yond using these two different boundary conditions. Finally 
we consider the effect of varying the magnetic diffusivity, 
rj. 

3. Results 

3.1. 2x2 case 

We take as our reference case the simulation starting from 
the checkerboard (2x2) initial condition, with vertical mag- 
netic field boundary condition at the upper boundary, hor- 
izontal resolution of 20.8 km, and rj = 11 km 2 s _1 . The 
magnetic Reynolds number of these simulations is below 
the critical value require d for small-scale dynamo action 
(|Vogler fc Schi issler 2007). This means that the simulations 
are too diffusive to show local dynamo action. Since our 
boundary conditions do not allow for incoming flux at the 
lower boundary, the magnetic field must inevitably decay. 
This can be seen qualitatively in Fig. [TJ which shows maps 
of the vertical magnetic field component on the horizontal 



surface z = 0, at the average (time and space) height of 
the optical (r = 1) surface. The temporal variation of the 
spatially averaged r = 1 height for this particular run was 
calculated for four different times and was found to have 
rms deviations of 8.5 km. For comparison the spatial rms 
deviations at any one time are about 4 times larger so that 
the r = 1 surface can be regarded as a corrugated suface 
whose average height only varies slightly. 





3w. 



Fig. 1. Time evolution of B z on the plane z = 0, corre- 
sponding to the average height of the optical surface (left 
panels) and the bolometric brightness (right panels) for the 
reference (2 x 2) case at 4 min (top), 80 min (middle), and 
200 min (bottom) after introduction of the initial field. The 
average unsigned magnetic field for these snapshots corre- 
sponds to about 200 G, 50 G, and 10 G, respectively. 



To quantitatively analyze the decay we consider the 
time evolution of the magnetic energy in a horizontal cut 
at a fixed geometrical height near z = 0, the evolution 
of the unsigned vertical magnetic flux through the same 
surface, and the magnetic energy in the whole computa- 
tional domain. We chose these three quantities because, 1) 
the vertical flux is the simplest quantity most relevant for 
comparison with observations, 2) the total magnetic energy 
because it gives information on the global properties of the 
simulation, and 3) the magnetic energy at the height of the 
vertical flux slice gives insight into the connection between 
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previous two. The evolution of these quantities is shown in 
Figure [U 
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Fig. 2. Decay of total magnetic energy in the entire simu- 
lation box, and of unsigned vertical magnetic flux at z = 
for the reference case. The thin lines are exponential fits. 
All curves have been normalized such that the fits have a 
value of unity at t = 0. 

An interesting feature of Figure [2] is that the unsigned 
flux, and the magnetic energy density at z = almost 
evolve in parallel. This runs against the naive expecta- 
tion that if the magnetic field decays at a rate a, so that 
B ~ ex.p(—at) then B 2 ~ exp(— 2at), suggesting that the 
magnetic energy should decay twice as fast as the unsigned 
magnetic flux. The reason why this relation does not hold 
can be seen qualitatively in Figure [TJ a significant amount 
of the magnetic flux is concentrated in a number of features 
and it is the number of such features, rather than the field 
strength in the features, which decreases in time. Therefore 
both the unsigned flux and the energy are proportional to 
the number of features and thus decay in parallel. This in- 
terpretation is confirmed in Figure [3j where we see that 
the percentage of surface magnetic energy corresponding 
to the different field strengths varies only weakly for the 
first 2 hours. In the third hour we begin to see an increas- 
ing contribution from weak field regions, as the number of 
strong- field elements diminishes. In particular we note that 
at 200 minutes there are only two reasonably strong mag- 
netic features, one of each polarity, left (cf. Fig. [TJ bottom 
panels) . 

Returning to Figure 2, the total magnetic energy in the 
entire computational domain decays faster than the sur- 
face energy densities and flux. We have estimated the de- 
cay rates by fitting an exponential function over the period 
from t = 30 min to t = 120 min. This period was chosen 
so that the initial phase, during which flux expulsion acts 
to increase the energy, is excluded. The e-folding time for 
the surface flux and magnetic energy, derived from the fits, 
is 70 minutes. The e-folding time for the magnetic energy 
over the entire domain is 52 minutes. The ratio between the 
decay rates of surface and volume integrated quantities is 
thus approximately 0.74, which is intermediate between the 
naive expectation of 0.5 and the expectation from magnetic 
elements of 1. This indicates that the magnetic structures 
become more volume filling (less intermittent) with depth. 
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Fig. 3. The percentage contribution to the total surface 
magnetic energy density per 100 Gauss, at different field 
strengths for different times (2x2 reference case). 



To characterize the decay in terms of a turbulent dif- 
fusivity, we consider the two-dimensional linear diffusion 
equation, dB z /dt = r] e tfV 2 B Zl where B z is the verti- 
cal component of the field on the plane z = 0. An as- 
sumption here is that the vertical flux diffuses as a pas- 
sive scalar. The slowest decaying eigenmode with the sym- 
metries of the initial condition has a spatial dependence 
sm(27rx/L) sin(27n//L) where L = GOOO/v^ km is the wave- 
length. The e-folding rates determined from the simulation 
can then be used to evaluate 7? e ff = 145 km 2 s _1 , which is 
about a factor of 13 larger than the explicit diffusivity in 
the simulation. 

The reconnection of the magnetic structures takes place 
when magnetic structures of opposite polarity come into 
close proximity. 

An example of such a reconnection event is shown in 
Figure |4] , where two opposite-polarity magnetic features 
have been brought into contact in an intergranular lane. 

The reconnection, which occurs above the surface, has 
locally heated the plasma. The strong down flows in the up- 
per photosphere are the result of the retraction of inverse- 
U loops (see Figure [5]) due to magnetic tension. Since the 
reconnection occurs close to the upper boundary, the corre- 
sponding U loops are rapidly removed from the system so 
that their signature in this reconnection event is very weak. 

A detailed description of the process invol ved and 
their o bservational signatures has been given by iDanilovid 
(|2009f ). The time scale for the individual reconnection 
events is of the order of 10 minutes, which is short com- 
pared to the decay rates found here. In what follows we 
concentrate on the statistical properties of the decay of the 
field rather than on the detailed reconnection processes. To 
do so we consider magnetic field lines, which are described 
by the function X, with 



dX(s,x ,yo,z ) _ B [X(s) 



ds 



B[X(s)]|+e 



(3) 



and 



X(s, xq, 2/ 0j zo) = (xo,y, 0, zo). (4) 
where e = 10 ~ 5 prevents problems where B = and 
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Fig. 4. Time series of a reconnection event (2x2 run, non-grey), approximately 30 minutes after introduction of magnetic 
field. The time between snapshots is 30 seconds. From top to bottom are the vertical component of the magnetic field, 
the vertical component of the velocity, the temperature and lastly the bolometric intensity. The first three quantities are 
shown at two heights, near z = km (lower rows) and near z = 300 km (upper rows). The horizontal extent of the maps 
is 1.75 x 1.75 Mm 2 . The intensity timeseries corresponds to the continuum intensity at 500 nm. The greyscale for the 
magnetic field covers the range from -530 G to +530 G at z = 300 km and from -1760 G to 1760 G at z=0; the velocity 
range is from -7 km/s to 7 km/s at z = 300 km and -9 km/s to 9 km/s at z = 0; the temperature range is from 4000° 
(blue) to 7500° (red) at z = 300 km and 5000° to 10000° at z — 0; and the intensity varies between 0.65 and 1.5 of the 
averge quiet Sun value. 



s is a coordinate along the field line, and xo,yo,zo is 
point in space corresponding to s = and the constant 
e = 10 _5 G prevents problems where B = 0. In regions 
where B ^> e, s is the length along the field line. 



For our analysis we consider the set of field lines pass- 
ing through the centre of each grid cell, i.e., we consider 
288 x 288 x 100 values of {xo, yo,zo}. This mean that near 
t = each cell will be sampled by 100 field lines, each 
of which will pass through the centre of at least one cell. 
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As time progresses stronger field regions will become even 
more oversampled and weak-field regions less sampled. At 
all times however each cell will be sampled by at least 
one field line. Through this we are able to capture cells 
which contain field lines entering and leaving the compu- 
tational box through different boundaries. For each field 
line, we integrate Eq. 3 numerically using a fixed step size 
As = 3.5 km until the field line exits the domain through 
either the upper or lower boundar^Q. Since the box is peri- 
odic in the horizontal directions almost all of the field lines 
pass through the upper or lower boundary. 




Fig. 5. The colored lines are a 3D visualization of field lines 
at beginning of the reconnection event shown in the SJ The 
grey-scale corresponds to the vertical component of the field 
at z = 0. The left panel is a side view, the right panel is a 
top view. 



The connectivity of the field line is defined by the 
boundary where it enters and leaves the computational box. 
A connectivity map is then created by considering the con- 
nectivity of all field lines. Since initially the field lines are 
all purely vertical, they can be classified as passing from 
the top to the bottom of the box or from the bottom to the 
top. As reconnection take place, some of these field lines 
are converted into U and inverse- U loops. Cells can con- 
tain field lines with different types of connectivities. 

The loops can then retract and leave through the bound- 
aries of the simulation domain. Figure [6] illustrates the con- 
nectivity of the field lines at three different heights for a 
time just after the initial period of flux expulsion and late 
in the simulation (at 16 and 142 minutes, respectively). 
The plane at height z = 420 km is in the upper photo- 
sphere, where the density has fallen considerably so that 
the magnetic tubes fan out and become smooth and space 
filling. The lowest layer, at z = —420 km, is within the 
convection zone, where the plasma motions dominate the 
magnetic field. Hence the magnetic field is more tangled 
and the connectivity map is more complicated. 

In the connectivity maps 16 minutes, the initial checker- 
board arrangement is still visible at all heights. Inverse- U 
loops, both ends of which connect through the lower bound- 
ary are present mainly near the borders of the opposing po- 
larities, indicating that reconnection has taken place. The 
U loops, both ends of which leave through the upper bound- 
ary mainly occur in the layers above the visible surface. The 
maps at 142 minutes are similar, except that the inverse- U 



1 The integration was limited to the range —3500km < s < 
3500 km. Since the computational domain is only 1.4 Mm in 
height most field lines left the box within this range. 




Fig. 6. Field line connectivities at 3 heights (from bottom 
to top: z = —420 km, z = and z = 420 km). As described 
in the text, field lines were traced from the centre of each 
of the 288 x 288 x 100 grid cells, so that on average each 
cell contains many field lines. The left panels correspond 
to t = 16 minutes, the right panels to t = 142 minutes of 
simulated time. White indicates areas where field lines pass 
from the bottom of the box to the top. Grey shows areas 
where field lines run from the top of the box to the bottom. 
Green indicates cells where all field lines have both legs 
passing through the bottom of the box (inverse- U loops). 
Blue indicates U loops with both legs passing through the 
upper surface. Light (dark) green indicates the presence 
of both inverse- U loops and field lines running from the 
bottom of the box to the top (top to bottom), within one 
grid cell. Light and dark blue regions are similarly defined 
but for U loops. The orange regions indicate where at least 
three types of field lines pass through a single grid cell. 



loops are present over an even wider region of space and 
the initial checkerboard pattern is less obvious. 

The predominance of inverse- U loops over U loops is a 
consequence of the fact that the reconnection occurs mainly 
above the optical surface, where the plasma beta is low 
and the magnetic field concentrations, expanding to become 
space filling, first come into contact with each other. The 
U loops then contain little mass and are quickly pulled 
upwards by the comparitatively strong tension force. The 
inverse-^/ loops, on the other hand, once they are pulled 
below the r = 1 level, are loaded with the mass from the 
subphotosphere and the tension force does not dominate 
over the other forces in this subsurface region. Magnetic 
tension is dominant in the current sheets in the low-/? region 
above the surface. 
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The retraction of the inverse- U loops is thus a much 
slower process. This explains the time evolution of the 
amount of flux in the different types of connectivities shown 
in Figure [71 Between t = 16 min and t = 142 min, the 
amount of unsigned flux corresponding to each type of con- 
nectivity at all heights has fallen substantially. The slowest 
decay is seen for the inverse- U loops, where the unsigned 
flux has fallen by a factor of approximately 3 below the 
r = 1 height and by a factor of 3-5 above. The unsigned 
flux associated with top-to-bottom and bottom-to-top field 
lines in the same time has fallen by a factor of approxi- 
mately 7, and the decrease in the flux of the U loops is 
even more dramatic. This reflects the fact that removal of U 
loops is a much faster process than the removal of inverse- U 
loops. Importantly, if the reconnection were mainly occur- 
ing in the sub-photospheric region, then the U loops would 
also be loaded with substantial mass and would only escape 
much more slowly, so that they would be more abundant. 




of opposite-polarity magnetic field. Again we quantify the 
decay rate by considering the decay of the total mag- 
netic energy in the computational domain. This is plot- 
ted in Figure [51 The e-folding time for the decay is found 
to be lhr 57 min. The dominant wavelength is 6000 km, 
so that the estimate for the turbulent diffusivity becomes 
77 eff 129 km 2 s _1 . 




2 3 
Time [hours] 

Fig. 8. Time evolution of the total magnetic energy starting 
from the 2x1 (striped) initial condition. The curve from the 
2x2 (reference) case is reproduced for comparison. 
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Fig. 7. Unsigned flux density as a function of height for the 
four connectivities. The solid lines correspond to the top- 
bottom and bottom-top connectivities, the dashed lines to 
inverse- U loops, and the dot-dashed line to U loops. The 
upper panel is for t = 16 minutes, soon after the initial tran- 
sient phase. Already a substantial amount of inverse- U loop 
flux has been generated. The lower panel is for t = 142 min- 
utes. The top-bottom and bottom-top connectivities use 
the same linestyle to reflect the B — ^ (— B) symmetry of 
the MHD equations. 



3.2. 2x1 initial condition. 

We now consider the value of rj e s based on the evolution 
from an initial condition consisting of two stripes (2x1) 



3.3. Upper boundary conditions 

We also consider the effect of varying the upper boundary 
condition on the decay rate. For the magnetic field we as- 
sume vertical or matched-to-potential boundary conditions; 
for the fl ows we consider either clo sed or open condit ions, 
following IsTein fc Nordlundl (|l998f ) and lBercikl (|2002[ ). We 
also vary the height at which the boundary condition is im- 
posed (z « 600 km above the r = 1 level to z « 840 km). 
In each case, we use the 2x2 initial condition as in the 
reference case. 

The decay of the total magnetic energy for the different 
cases can be seen in in Figure [9] The e-folding time for the 
magnetic energy and the estimated turbulent diffusivities 
are listed in Table 1. The values for the turbulent diffu- 
sivity resulting from the use of different upper boundary 
conditions vary by a factor of 3, from 145 km 2 s _1 for the 
vertical field, closed boundary applied at z « 600 km to 
341 km 2 s _1 for the potential field, open boundary applied 
at z w 840 km. 



3.4. Explicit magnetic diffusivity parameter, rj 

We investigated the role of the explicit value of the in- 
put parameter rj by performing a simulation with twice 
as many grid points in the horizontal directions and with 
77 = 5.6 km 2 s _1 . The initial condition was the same as 
in the reference case, and the time evolution is shown in 
Figure [TUJ The fit to the total energy corresponds to an e- 
folding time of lhr 15 min, and hence an effective diffusivity 
?7eff = 101 km 2 s _1 . We note that the the fit to the curve in 
Figure [10] is sensitive to the time interval over which the fit 
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Table 1. Summary of Magnetic energy decay rates and r] e ^ 
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Fig. 9. Time evolution of the total magnetic energy for the 
2x2 case with several different upper boundary conditions 
located at one of two heights. The 'normal, closed, ver- 
tical' case refers to our standard 2x2 run. The 'normal, 
closed, potential' refers to a case where the magnetic field 
was matched to a potential field at the top boundary. For 
the 'tall' cases the top boundary was moved upwards by 
280 km, and 'open' refers to cases where flow, tempera- 
ture and pressure were smoothly extrapolated through the 
upper boundary. 

was made. Fitting over the interval 30-180 minutes leads 
to essentially the same estimate for ?7 e ff as was obtained for 
the 2x2 reference case. In this regard it is worth noting 
that the late-stage evolution in most cases is different from 
that in the first 2 hours or so. This is presumably due to 
the fact that the random motions move energy from those 
present in the initial condition modes (e.g. 2x2) into other 
modes (e.g. 1x2 or 3x2). When these modes dominate the 
decay, a different decay rate is to be expected. 

4. Discussion 

The values of r] e ^ derived from the various runs and listed 
in Table 1 are within the range deduced from observations 
over a much larger spatial e xtent and longer time period 
(see lSchrijver fc Zwaanll2000[ Table 6.2). To understand the 
differences between the estimated values of 77 e ff (i.e. the de- 
cay rates) for the different simulation setups, we note that 
the removal of flux in these simulations is primarily deter- 
mined by the rate at which opposite polarity patches are 
brought close enough to one another for the explicit dif- 
fusivity to become effective. The differences between the 
values of 7? e ff can be understood in terms of how close the 
magnetic polarities need to come for rapid diffusion to set 



Fig. 10. Time evolution of the total magnetic energy for 
a checkerboard (2x2) run with r] = 5.6 as a function of 
time. The same curve from the reference case is shown for 
comparison. 



in. In loose terms this can be rephrased in terms of a 'col- 
lisional cross section' of the magnetic elements. When we 
use the potential boundary condition, the magnetic field 
expands rapidly near the upper boundary, and hence mag- 
netic elements which are still separated near the surface 
can still come into contact in the upper layers. This is the 
reason why the reconnection mostly takes place above the 
surface. It also explains, qualitatively, why the potential 
boundary condition leads to a higher r] e ff then was found 
in the reference case. 

We may compare the values of r] e ^ with the estimate, 
^turb ~ ^rms^corr/3, for the diffusion of a passive scalar 
subject to random motions. In this estimate v rms is the rms 
velocity of the random motions and L corr is the correlation 
length scale of the motions. The first question is then which 
motions to consider. Since we are mainly concerned with the 
horizontal motions of the magnetic flux tubes we consider 
only the horizontal velocity field in a z = slice from the 
simulations. The rms horizontal velocity is then 2.4 km s _1 . 
However, for the horizontal transport of the flux concen- 
trations in the intergranular lanes, the random horizontal 
motion of the granules is the relevant quantity. The corre- 
sponding rms value, for example as determined by local cor- 
relation tracki ng, is about a factor of two smaller than the 
above value dRieutord et al.ll200ll : iGeorgobiani et al.ll2007t 
Matl och et al. ,1 120101) . ie Vr™* = 1.2 km s 1 . For the correla- 
tion length of the horizontal velocity, we take L corr = 1 Mm, 
which is about the size of a granule. These values yield 
?7 e ff = 400 km 2 s _1 . This very simple, approximate value 
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is slightly above the actual range of values found in the 
simulations. 

5. Conclusions 

We have studied the decay of simple arrangements of mag- 
netic field in the near-surface layers of the Sun. We found 
that magnetic field is removed from the surface when mag- 
netic elements of opposite polarity, advected by the gran- 
ular flows, come into close proximity above the surface. 
Reconnection then produces U loops, which quickly escape 
through the upper boundary, and inverse- U loops, which 
retract more slowly into the subphotosphere. 

During most of the decay phase, the relative contri- 
butions of magnetic energy coming from different field 
strengths are largely independent of the amount of flux. A 
weak field tail, however, does appear late in the evolution, 
when the unsigned flux levels are very low. 

The decay of the magnetic field is reasonably well de- 
scribed by turbulent diffusion concept with plausible values 
of Tfeff . The value of rj e fi was found to depend on the upper 
boundary condition. The potential-field, open-boundary 
condition applied to the tall box is probably the most real- 
istic in this regard. This case has r^ff ~ 340 km 2 s _1 . The 
range of values from the different simulations, from 100 to 
340 km 2 s _1 , gives an indication of the uncertainty due to, 
for example, the choice of the upper boundary condition. 
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